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Abstract. In the field of Energy Density Functionals (EDF) used in nuclear structure 
and dynamics, one of the unsolved issues is the stability of the functional. Numerical 
issues aside, some EDFs are unstable with respect to particular perturbations of the 
nuclear ground-state density. The aim of this contribution is to raise questions about 
the origin and nature of these instabilities, the techniques used to diagnose and prevent 
them, and the domain of density functions in which one should expect a nuclear EDF 
to be stable. 
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The Energy Density Functional (EDF) method is a broadly used tool in the study 
of nuclear structure and dynamics, allowing to describe medium-mass to heavy nuclei 
with reduced computational cost. This is achieved thanks to the implicit resummation 
of quantum many-body correlations into an effective energy functional of local densities, 
which is, in principle, minimized for the intrinsic density of the many-body ground state 
(i.e. the density in the center-of mass frame PQ, possibly breaking several symmetries of 
the Hamiltonian) . Non-relativistic EDFs mostly used in nuclear physics until now are 
of the Skyrme type, i.e. semi-local functionals with a postulated analytical form and no 
clear connection to the "bare" two- and three-nucleon interactions. 

The Skyrme EDF has been employed not only in static "mean-field" calculations, 
but also for the study of small- or large-amplitude collective dynamics in schemes 
such as Time-Dependent-Hartree-Fock (TDHF), or Quasiparticle Random-Phase 
Approximation (QRPA). Furthermore, the power of the Multi- Reference EDF method 
[21 [3], which includes collective- motion correlations into ground and excited states 
through the resolution of the Hill- Wheeler equations, has been thoroughly demonstrated. 

In the present article we concentrate on "Skyrme-like" or semi-local functionals. 
One of the intrinsic, unsolved issues of these functionals is the stability and convergence 
of calculations performed with them. Besides numerical issues linked with finding an 
energy minimum, especially a constrained one, some functionals appear unstable with 
respect to specific perturbations of the nuclear density. This instability means that self- 
consistent calculations will yield unphysical configurations at energy minimum, break 
symmetries expected in the minimum-energy density (i.e. other than those supposed to 
be explicitly broken such as spherical symmetry) and/or yield a divergent energy. This 
is now understood to be an intrinsic defect of some functionals, independently from 
numerical details [HE]. 

The Skyrme EDF can be expressed through the local energy density %, 

H = £-T + H + Hi, (1) 
2m 

where isoscalar (t = 0) and isovector (t — 1) terms have been separated; both have the 
same structure, i.e. for even-even nuclei, 

H t = C p tP 2 t + Cl Pt T t + C^pApt + C7 J Pt VJ t + C J t J 2 t , (2) 

where p, r, and J are respectively the matter, kinetic, and spin-current densities, 
an index t = 0, 1 indicating respectively isoscalar and isovector quantities; standard 
definitions of the latter can be found in Refs. [21 E]- in typical Skyrme parametrizations, 
density dependence is introduced as Cf = C^+C^Pq, with p = p n +p p . Generalizations 
involving a dependence of coupling constants on the isovector density p\ = p n — p p 
and/or more involved functions of densities are being studied. Recent work [7J [8] 
aiming at building non-empirical, yet semi-local, functionals through the Density Matrix 
Expansion (DME) method [9] also leads to such generalized forms. 

The instability of Infinite Nuclear Matter (INM) and nuclei with respect to global 
spin and/or isospin polarization is a common problem for functionals derived strictly 
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Figure 1. Radial density profile in nuclei 40 Ca (top panel) and 56 Ni (bottom 
panel) obtained through EDF calculations in spherical symmetry, respectively with 
functionals LNS [M] and SkP [13] . Each curve corresponds to the number of iterations 
indicated in the legend. Stable results obtained with SLy5 [15] are given for reference. 



as the expectation value of a Skyrme effective interaction, and has been repeatedly 
discussed and addressed [TQl [TTJ [12] . However, these pathologies of the INM Equation 
of State (EOS), usually occurring at high densities, hardly affect nuclear structure 
calculations in practice. 

More recently, a divergence of the total energy and densities has been observed 
in EDF calculations of finite nuclei [4J performed with specific parameterizations 
of the Skyrme functional [T3l [T3] . Fig. [1] displays densities obtained with these 
parameterizations in two sample nuclei at various (large) numbers of iterations of the 
self-consistent equations. Characteristic traits of this phenomenon include unnatural 
configurations (large- amplitude spatial oscillations of densities, with separation of 
protons and neutrons) and unnaturally large (negative) total energy. This instability 
was analyzed by studying the reponse of INM, used as a model system, to finite-size 
density perturbations at the Random Phase Approximation (RPA) level. 

When calculated with the central terms of the residual interaction deduced from 
a semi-local functional [16], the analytic expressions for the RPA response functions 
are rather involved, even more so when the spin-orbit [17] and tensor [18] terms are 
added. For the sake of discussion, it is useful to consider the case of a functional with 
Cl = Cf = C t VJ = (which is strictly local, despite the remaining pAp, "finite-range" 
terms) considering isospin density fluctuations only. One then finds the textbook result 
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for the response function in symmetric INM, 
jrf \ 4n (w,q) 

H[u},q) = a 3 

l-8[Cf-C 1 A V]n (o;,q) 

where u is the excitation energy, q the transferred momentum, or wave number of 

the density fluctuation, H(u>, q) the response function, and U (cu, q) the non-interacting 

response, or Lindhard function [19]. The value H(lo = 0,q) corresponds to the static 

susceptibility of the system to finite-size perturbations. It should be positive (with 

the above sign convention) for all values of q and the density p$. (We recall that the 

coupling constants in the above equation can depend on the latter, as does n .) A 

change of sign with either variable is accompanied by a pole signaling the existence of 

a zero-energy collective mode. We see that the short-wavelength (high q) behaviour is 

driven by the coefficient C^ p , whose value correlates with the occurrence rate of the 

instability observed in calculations of finite nuclei [4]. A similar convergence issue was 

identified in Ref. [20], and confirmed in Ref. [5], concerning the susceptibility of the 

functional SkO [21 J to fluctuations of the spin density. 

We could categorize these instabilities in terms of the importance and role of the 

finite size of the system under consideration. First of all, Landau's theory of Fermi 

fluids [22] gives quantitative constraints for the stability vs. long-range perturbations 

of INM in terms of Landau parameters. Second, the more general random-phase 

approximation (RPA) calculations in INM can, at least qualitatively, link zero-energy 

modes to oscillations of densities in static calculations of nuclei. Third, in finite nuclei, 

instabilities originating at the surface may not be amenable to a diagnosis through INM 

calculations. This may be emphasised, in particular, by the recent attention brought to 

more involved density dependences [23]. This leads to ask a first open question: 

1 Can we quantitatively understand the link between instabilities in INM and finite 
nuclei? In particular, what INM densities are relevant to the study of nuclei? How 
does finite size modify the results from INM? 

The principal manifestation of instabilities described up to now seems to be 
divergence of given quantities. However, let us point out that this is not mandatory 
for a functional to exhibit unstable behaviour. A schematic representation of the EDF 
energy with respect to an abstract collective coordinate x in the EDF calculation (e.g. in 
Fig. [I} the magnitude of the density oscillations) is given in Fig. |2j Here the dot at x = 
represents the physically expected configuration. Curve (a) corresponds to the (desired) 
situation with a stable configuration, whereas curve (d) represents a well-developed 
instability. Curves (b) and (c) represent situations with a metastable minimum which 
can be made to evolve to more energetically favored configurations by some perturbation. 
In case (b) a true, global variational solution exists; whereas in case (c) it does not. Here, 
RPA (i.e. linear response) analysis would only warn about situation (d). 

Case (b) was encountered in Appendix B of Ref. [21], where a parameterization with 
large negative values of C/ was found to yield, in seemingly random nuclei, an unphysical 
shell structure. The plot of the energy landscape as a function of the magnitude of the 
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Figure 2. Schematic picture of energy landscape with respect to collective coordinate 
x, assumed to be constrained while other degrees of freedom are variationally 
optimized, for stable (a) and various unstable cases. See text. 



spin-current density J in Fig. 35 of Ref. [21] indeed corresponds quite well to case 
(b), with two minima, the unphysical one being energetically favored. Thus, it seems 
appropriate to ask: 

2 How can we diagnose instabilities not accessible with a small-amplitude expansion 
around the expected ground state (as in RPA)? 

In Ref. [25]; a TDHF calculation of a central 16 0- 16 reaction done with the 
functional SkM* [26] showed that after contact, the two nuclei separated exhibiting a 
toroidal spin polarization. This puzzling result might be analyzed from the standpoint of 
stability analysis: is this configuration lying lower in energy than a more "reasonable" 
one? If so, this result might well have to be discarded. It is also possible that a 
treatment beyond the static considerations above (i.e. referring to energy minimization) 
be necessary for a thorough understanding, and hence: 

3 What are the manifestations of instabilities developing on top of excited states, 
probed in dynamical calculations? Can we reduce their study to a model problem, 
such as the INM analysis above? 

Finally, current efforts tend towards producing non-empirical semi-local functionals 
through the DME applied to low- momentum interactions [2Z] ■ The DME, more than a 
simple range expansion, consists of a series of finite-size corrections to the INM density 
matrix [9j [7J [8]. It can thus be expected to be accurate for smoothly varying densities 
such as those occurring naturally in nuclei. However, there is no guarantee that this 
accuracy will be preserved for pathologic configurations such as in Fig. [TJ In other 
words, extra care should be taken to verify that the DME functional does not exhibit 
pathological behaviour beyond its expected domain of validity. Indeed, by going from 
an exact Hartree-Fock functional deduced from a realistic Hamiltonian to a semi-local 
functional, spurious high-momentum components are introduced e.g. in the residual 
interaction and response functions (Eq. (3)). The apparent link between RPA linear 
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response analysis and self-consistent EDF calculations suggests that the "singularity" 
due to the semi-local nature of the functional can make divergences, naturally expected 
in calculations beyond first order of perturbation theory performed with a zero-range 
interaction, creep into "mean-field" ones through self-consistency. Our final question is 
thus: 

4 How far in density should we expect the next-generation nuclear EDF to be stable, 
and in what domain of variation for the reference Slater determinant and associated 
densities? In other words, shouldn't semi-local functionals be specified with a well- 
defined "cutoff" or resolution scale? 

It is especially important to understand this issue for functionals containing higher- 
order derivative terms (yielding higher orders of q in Eq. (3)) such as proposed 
in Ref. [28], whether they are built from an underlying Hamiltonian or purely 
phenomenological. 

In summary, despite their great usefulness, semi-local nuclear energy density 
functionals of the Skyrme type come with several kinds of stability-related issues, which 
can compromise the significance of results obtained with them. Indeed, whereas in 
the process of developing new nuclear functionals, one naturally focuses on minimum 
(ground-state) energies and associated density configurations, these functionals will be 
probed in different ways in actual applications, be it during energy minimization or in 
dynamical extensions of the theory. These issues should thus be fully understood and 
kept in mind for future development. 
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